For weather data…
For health outcomes data, we are interested in the incidence of asthma, skin cancer and lung cancer across time and space. We obtained age-adjusted incidence rates data for these three health outcomes for the U.S. states New York, Ohio and Pennsylvania over multiple years. The following bullet points provide links from where we obtained the data.
We also obtained age-adjusted incidence rates data for these three health outcomes for the U.S. states New York, Ohio and Pennsylvania at the county-level at fixed time points. The following bullet points provide links from where we obtained the data. The asthma data is provided for the year 2016 and the skin cancer and lung cancer data is averaged over the years 2014-2018.
Each of these data sets were exported and each imported into R. For each of the three health outcomes, the data sets were merged across the three states for the longitudinal data and for the cross-sectional data. This effort resulted in six data sets. The data import and cleaning of the original data can be found in the data import folder on our github repository.
We will now import the six data sets and merge all longitudinal data together and merge all cross-sectional data together to have two data sets.
First, let’s combine the longitudinal lung and skin cancer data at the state level.
# read in longitudinal state data for lung and skin cancer
# longitudinal state level lung cancer data
state_lc_data = read_csv("./data/lung/state_lc_data.csv")
# longitudinal state level melanoma data
state_mel_data = read_csv("./data/incidence_melanoma_data/state_melanoma_age_adjusted.csv") %>%
rename(c("outcome" = "cancer_type"))
# combine cancer outcome data
lc_mel_state = bind_rows(state_lc_data, state_mel_data) %>%
select(state, year, outcome, age_adjusted_incidence_rate)
Next, let’s combine the longitudinal asthma data with the longitudinal cancer data at the state level.
We begin by reading in the longitudinal state-level asthma data.
# Read in longitudinal state-level asthma data
# Rename and select columns
asthma_state <- read_csv(here::here("data", "asthma_state.csv")) %>%
rename(age_adjusted_incidence_rate = rate) %>%
select(state, year, outcome, age_adjusted_incidence_rate)
Next, we can combine the longitudinal asthma and the longitudinal cancer state-level data.
# Bind rows of longitudinal asthma and longitudinal cancer by state and arrange
# accordingly
lc_mel_asthma_state <- bind_rows(lc_mel_state, asthma_state) %>%
arrange(state, year, outcome)
Lastly, let’s write out the combined data.
write_csv(lc_mel_asthma_state, here::here("data", "lc_mel_asthma_state.csv"))
First, let’s combine the cross-sectional lung and skin cancer data at the county level. We will add a column here for the county FIPS code for mapping purposes.
# read in cross-sectional county data for lung and skin cancer
county_lc = read_csv("./data/lung/county_lc.csv")
county_mel = read_csv("./data/incidence_melanoma_data/state_county_melanoma_incidence_2014_2018.csv")
fips_codes = read_csv("./data/fips_codes.csv")
# cleaning to make the data sets compatible
county_mel = county_mel %>%
rename(c("outcome" = "cancer_type")) %>%
drop_na()
county_mel <- county_mel[!(county_mel$county == "Ohio" | county_mel$county == "Pennsylvania" |
county_mel$county == "New York"), ] %>%
filter(!grepl('SEER', county))
county_mel$county <- gsub(" County","", county_mel$county)
county_mel$county[county_mel$county == "St Lawrence"] <- "St. Lawrence"
county_lc$county <- gsub(" County","", county_lc$county)
# binding rows
lc_mel_county = bind_rows(county_lc, county_mel) %>%
filter(complete.cases(.))
# Add FIPS column to the data
fips = fips_codes %>%
janitor::clean_names() %>%
rename(c("county" = "name"))
lc_mel_county = merge(lc_mel_county, fips, by = c("state", "county"))
Next, we will read in the cross-sectional asthma county data and combine it with the cross-sectional cancer data.
We begin by reading in the cross-sectional asthma county data, renaming column names and adding a column for county FIPS code.
# Read in cross-sectional asthma data
# Only get county data for the year 2016 since we need county data at
# a fixed time point (i.e. not longitudinal)
# Adjust column names and types in prep for merge
asthma_county <- read_csv(here::here("data", "asthma_county.csv")) %>%
filter(year == 2016) %>%
select(-year) %>%
rename(age_adjusted_incidence_rate = rate) %>%
mutate(age_adjusted_incidence_rate = as.double(age_adjusted_incidence_rate))
asthma_county <- merge(asthma_county, fips, by = c("state", "county"))
Next, let’s combine the cross-sectional asthma and cancer county data.
# Bind rows of cancer and asthma by county and arrange accordingly
lc_mel_asthma_county <- bind_rows(lc_mel_county, asthma_county) %>%
arrange(state, county, outcome)
Lastly, let’s write out the combined data.
write_csv(lc_mel_asthma_county, here::here("data", "lc_mel_asthma_county.csv"))
Now, that we have combined the health outcomes data, we can explore.
To discuss how climate change has impacted the rate of certain negative health outcomes, let’s first assess their general prevalence.
Clearly, according to the map below, lung cancer is a relevant issue in our three states of interest (NY, PA, and OH). They also provide a good range of examples, as we can see multiple counties with high age-adjusted incidence rates in Ohio, a few in New York, and virtually none in Pennsylvania.
# plotly of county level data
county_map_lc
We can also check how these rates have changed by state over time:
# plotly of county level data
fig
state_mel_plot

Let’s plot the combined longitudinal state-level health outcomes data with year on the x-axis, age-adjusted incidence rate on the y-axis, stratified by state (color) and outcome (facet). Please note each y-axis has its own scale.
# Scales = free allows each y-axis to be independent of the others
lc_mel_asthma_state %>%
ggplot(aes(x = year, y = age_adjusted_incidence_rate,
color = state)) +
geom_line() +
facet_wrap(~outcome, scales = "free") +
theme_bw() +
labs(
title = "Age-Adjusted Incidence Rate for each Outcome by State",
x = "Year",
y = "Age-Adjusted Incidence Rate",
color = "State"
) +
theme(plot.title = element_text(hjust = 0.5))
